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Influence of defects on ferroelectric and electrocaloric properties of BaTiOa 
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We report modifications of the ferroelectric and electrocaloric properties of BaTiOs by defects. 
For this purpose, we have combined ab initio-hased molecular dynamics simulations with a simple 
model for defects. We find that different kinds of defects modify the ferroelectric transition tem¬ 
peratures and polarization, reduce the thermal hysteresis of the transition and are no obstacle for 
a large caloric response. For a locally reduced polarization the ferroelectric transition temperature 
and the adiabatic response are slightly reduced. For polar defects an intriguing picture emerges. 
The transition temperature is increased by polar defects and the temperature range of the large 
caloric response is broadened. Even more remarkable, we find an inverse caloric effect in a broad 
temperature range. 

PACS numbers: 


1. INTRODUCTION 

As conventional cooling is presently one of the main 
carbon dioxide sources, more advanced technologies are 
urgently needed. In the last years the electrocaloric effect 
(ECE) has come into focus as a promising new cooling 
mechanism.Et^ Although the change of temperature in a 
ferroelectric material due to a varying electrical field is 
known for decades, giant caloric responses up to 12 K 
have been obtained only recently.^] The understanding of 
the ECE is still quite unsatisfactory. Materials with a 
large and reversible adiabatic response in a proper tem¬ 
perature range, and which are ecologically save and eco¬ 
nomically viable, have to be found. 

So far little is known about possible limitations of the 
ECE due to defects. However, defects such as missing 
oxygen atoms or impurities and dopants are common 
in ferroelectric materials. Even though a large ECE of 
0.5 K/300 kV/cm has been found for doped BTO in an 
industrially manufactured multilayer capacitor,!^ we are 
not aware of any systematic investigation on the influence 
of doping and defects on the ECE. The present paper 
aims to fill the gap in the understanding and optimiza¬ 
tion of the ECE in the presence of defects. 

Different experimental and theoretical studies deal 
with the modification of ferroelectrics by doping and 
defects.l^^^Both, doping and defects, influence remanent 
polarization and polarization switching, which are impor¬ 
tant for the caloric response. Evidently, even the modifi¬ 
cation of the ferroelectric phase diagram under such dop¬ 
ing is still unclear. For example, reduction, insensitivity, 
and increase of the ferroelectric transition temperature 
in BaTiOs (BTO) due to doping have been found.l^^tlil 
In the present paper we use molecular dynamics simula¬ 
tions employing an effective Hamiltonian based on ab ini¬ 
tio parametrization to improve the understanding of the 
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ferroelectric phase diagram by defects. The same model 
has been used successfully in order to mod el ferroele ctric 
phase diagrams and the ECE in literature.l^®^E^E^ 

In particular, experimental and theoretical investiga¬ 
tions reveal the important role of transition me tal dopin g 
on fatigue behavior in ferroelectric perovskites. b^b jt 
has been found that an internal bias field builds up in 
field polarized (poled) Cr doped BTO, which could be 
related to the alignment of defect dipoles with the over¬ 
all polarization.li2 in this context Erhart and co-workers 
performed detailed ab initio based transition state theory 
simulations on transition metal doped perovskites.ff^ff^ 
They found that polar 0-vacancy dopant complexes, 
which induce local dipoles, form immediately. It is most 
favorable if these dipoles align parallel to the overall po¬ 
larization. However, their relaxation is rather slow and 
the defect dipoles thus give rise to aging and fatigue, i.e. 
broadening and shift of the held hysteresis. 

Also, these changes of the phase diagram may largely 
modify the caloric response and its reversibility. For the 
first time, the change of operation range and magnitude 
of the ECE with different kinds of defects is investigated 
in the present work. Also the influence of different mea¬ 
suring protocols is revealed. With measuring protocol 
we mean the influence of the sample history, i.e. poled 
or as-prepared samples, cooling-down or heating-up sim¬ 
ulations and the use of unipolar and bipolar electrical 
fields. 

We show that a reduced local polarization reduces Tc, 
polarization, and the ECE, while local defect dipoles with 
slow relaxation dynamics and large dipole moment in¬ 
crease Tc and may even induce an inverse caloric re¬ 
sponse. The paper is organized as follows. First, the 
computational methods and the model for defects are 
discussed in Sec. [21 A more detailed discussion on the 
accuracy of the method and a summary of convergency 
tests can be found in Sec. jH In Section [3] the influence 
of non-polar defects on the ferroelectric phase diagram 
and the ECE are discussed. Section [i] deals with the in¬ 
fluence of polar complex dipoles on the phase diagram of 
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FIG. 1: (Color online) Schematic sketch of the used defect 
model, (a) non-polar defects with locally frozen soft mode 
amplitude (b) fixed defect dipoles to which parallel or anti¬ 
parallel external electric field is applied. 


BTO and the ECE for different defects and measuring 
protocols. Finally, conclusions and outlook are given in 
Sec. [Sl 


2. COMPUTATIONAL MODEL 

Molecular dynamics (MD) simulations have been 
performed employing the feram code http://loto. 
sourceforge.net/ferEmi/ developed by Nishimatsu et 
al^ based on an effective Hamiltonian 

J[/f* _ A//* _ 

^ ^ . 2 ^ ^ 

R,a 

+ ryi, •••,%) + {m}), 

( 1 ) 

with rji,... ,r]Q the six components of homogeneous strain 
in Voigt notation. is the self energy of the 

local mode, is the long-ranged dipole-dipole 

interaction, is the short-ranged interaction 

between local soft modes, ..., r^g) is the 

elastic energy from homogeneous strains, ‘''''“({tu}) 
is the elastic energy from inhomogeneous strains, 
^coup.homo|'|.j^|^ ^ ^ is the coupling between the 
local soft modes and the homogeneous strain, and 
in*'°({M}, {ic}) is the coupling between the soft 
modes and the inhomogeneous strains. 

Instead of treating all atomic positions as degrees of 
freedom, the collective atomic displacements are coarse¬ 
grained by local soft mode vectors u{R) and local acous¬ 
tic displacement vectors w{R) of each unit cell at R in 
a simulation supercell. Therefore, the number of de¬ 
grees of freedom per unit cell is reduced in a first step 
from 5 atoms times the three cartesian directions (= 15) 
to 2 three-dimensional vectors, .^dipoie/^ Sr a 
and K{R) are the kinetic energies pos¬ 

sessed by the local soft modes and the local acoustic dis¬ 
placement vectors along with their effective masses of 
M^jpoie and . In a second coarse graining step 

an internal optimization of E{w{R)) is used reducing 


the degree of freedom to the three components of the 
soft mode u{R). Details of this Hamiltonian are ex¬ 
plained in Refs. The set of parameters for the 

effective Hamiltonian for BTO have been obtained by 
density functional theory simulation at T = 0 K and are 
listed in Ref. |23l 

Within the molecular dynamics (MD) simulations, pe¬ 
riodic boundary conditions and a cell size of 96 x 96 x 96 
are used unlike otherwise stated. This corresponds to 
884736 f.u. i.e. 4423680 atoms. For cooling-down and 
heating-up simulations on a dense temperature grid of at 
most 5 K steps, the pre-converged configuration of previ¬ 
ous steps is used as input for each temperature step and 
thus the equilibration time can be reduced to 40.000 fs. 

The procedure of direct simulation of the ECE involves 
three steps: first, we equilibrate the system in an external 
electrical field, E along [001], constrained to a constant 
number of particles (N), constant pressure (P), and con¬ 
stant temperature (T) (NPT). For this purpose we apply 
the Nose-Poincare thermostat.l^Next, we change to adi¬ 
abatic conditions keeping the number of particles N, the 
pressure (P), and the total energy (E) constant (micro 
canonical NPE ensemble) and let the ensemble evolve in 
time by the leapfrog method. Simultaneously, the field 
is ramped down with a rate of 0.05 kV/cm/fs, see Sec.[^ 
for a detailed discussion. The final state at the end of 
the constant temperature MD is used as the initial state 
for the constant energy MD. As last step we monitor the 
kinetic energy after 10.000 further equilibration steps. A 
time step of at most At = 2 fs is used in both ensembles. 

In the Hamiltonian, only the three components of the 
soft mode vector are explicitly taken into account. Thus, 
the number of degrees of freedom is reduced from 15 to 
3 which in turn reduces the specific heat. As a result the 
adiabatic response obtained within our model is overes¬ 
timating the response of real BTO. However this overes¬ 
timation is an intrinsic feature of the model and thus the 
qualitative trends for different kinds of defects are not 
affected. In addition, it has been shown that the leading 
error of this deviation ca n be c orrected with a rescaling of 
the temperature by 15/3,1^^^ see Sec. 6. In the following 
all results are rescaled correspondingly in order to get an 
estimate of the actual temperature change. 

Defects are modeled by freezing the local soft mode on 
random positions to a certain value. A local reduction of 
polarization by a defect is modeled for the limiting case of 
non-polar defects. Therefore, the local soft mode ampli¬ 
tude in some cells is frozen to zero, see Fig. [^(a). Defect 
dipoles with slow relaxation are modeled by dipoles with 
fixed amplitude and direction, see Fig. [^(b). We assume 
a ferroelectric coupling between defect dipoles and the 
surrounding BTO matrix. Thus, we investigate the in¬ 
fluence of the modified dipole-dipole interaction by local 
changes of the polarization. The used coarse graining 
procedure does not allow to model atomic relaxations 
near the defect. In addition, we neglect the direct cou¬ 
pling of defects to the external electrical field and pertur¬ 
bations of the ferroelectric matrix near the defect, e.g.. 
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FIG. 2: (Color online) Polarization as function of tempera¬ 
ture of BaTiOs; Red: heating-up simulation; Blue: cooling- 
down simulation; Black: the sample has been equilibrated 
in 100 kV/cm and polarization and temperature have been 
recorded after field removal (ECE protocol). In all cases the 
polarization along [100],[010] and ]001] is shown. 

For the ECE protocol the data without pressure correction 
(dashed black lines) is opposed to results for an external pres¬ 
sure correction of p = —O.OOST, treating the acoustic degrees 
of freedom explicitly, and shifted by 150 K to lower tempera¬ 
tures (black stars). 


by a modified volume or modified short range interac¬ 
tions. 

For further details on the method and its excellent 
reliabilit y for the mod eling of the ferroelectric phase 
diagramj ^^ l ^^ l the or the influence of small 

distortions to the ideal syste m e .g., due to strain, 
alloying,!^ or finite-size effectd^^lMI refer to literature. 


3. BATIO 3 AND NON-POLAR DEFECTS 

The simulated polarization as function of tempera¬ 
ture of BTO is illustrated in Eig. At high tem¬ 
peratures BTO is in the paraelectric cubic perovskite 
structure. Under cooling, three coupled structural- 
ferroelectric phase transitions exist, to the tetragonal, 
orthorhombic, and rhombic phases, with spontaneous po¬ 
larization along [100], [110], and [111], respectively. We 
find transition temperatures of 275 K, 150 K, and 70 K 
for cooling-down simulations in qualitative agreement 
with the experimentally found phase sequence. In ad¬ 
dition, spontaneous polarization and strain of all phases 
as well as the first order character of all transitions are 
as in experiment.^ The polarization jumps by about 
30 /iC/cm^ at Tc and a pronounced thermal hysteresis is 
present between cooling-down and heating-up because we 
neglect nucleation sites such as defects, surfaces, or grain 
boundaries, see Eig. Experimentally, a thermal hys¬ 
teresis of about 10 K has been found for BTO ceramics.l^ 
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FIG. 3: (Color online) Polarization of BaTiOs as function 
of temperature T with and without non-polar defects. Defect 
concentrations ranges from 0-3%. Upward and downward ar¬ 
rows indicate cooling-down (blue) and heating-up (red) sim¬ 
ulations, respectively. Polarization under an external held of 
100 kV/cm along z is given in black. As the external held is 
larger than the critical held strength, the held induced tran¬ 
sition does not show thermal hysteresis anymore. 

In addition to cooling-down and heating-up simulations 
the polarization profile as obtained in the ECE measure¬ 
ment is included in Eig.j^in black. Eirst, the material has 
been equilibrated in an external field of 100 kV/cm, and 
after the field has been removed, polarization and tem¬ 
perature have been sampled. The polarization profile in 
this case mainly coincides with the heating-up simula¬ 
tions. 

Despite this excellent qualitative description of the 
material, all transition temperatures are underestimated 
by about 100 K within the model and parametrization 
used in our study, cf. Refs. I22I23I This may mostly 
by attributed to the fact that thermal expansion is 
strongly underestimated in the model. One can partly 
compensate for this with a semi-empirical effective 
negative pressure given by p = —0.005T GPa.l^ By this 
approach, Tc of the paraelectric to ferroelectric transi¬ 
tion increases to about 411/436 K under cooling-down 
and heating-up simulations, thus slightly overestimating 
experimental values. Apart from this temperature 
shift, the obtained phase diagram is hardly modified. 
Especially, the polarization profile as obtained under the 
ECE protocol without pressure correction (black dashed 
line) and the data with pressure correction shifted by 
150 K ( black stars) show qualitatively the same trend, 
see Eig. [2p^ All results presented in the following have 
been obtained without the semi-empirical correction 
unlike otherwise stated. However, all qualitative trends 
have been cross-checked with the effective pressure 
correction. 

The ECE in BTO can be understood by simple ar¬ 
guments. Dipole ordering and magnitude of the local 
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FIG. 4: (Color online) Change of polarization for different 
concentrations of non-polar defects after an external field of 
100 kV/cm has been removed. 
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FIG. 5: (Color online) Approximate adiabatic temperatnre 
change if an external held of 100 kV/cm is removed from sam¬ 
ples with different concentration of non-polar defects. Dotted 
lines illustrate the linear reduction of peak position and mag¬ 
nitude of the ECE. 


polar off-centering increase if an electrical field is ap¬ 
plied, see Fig. The induced polarization vanishes, af¬ 
ter the field is removed, see Fig. This reduction is 
related to an increase of the disorder of the local dipoles, 
i.e. an increase of the entropy under isothermal condi¬ 
tions, cf. Fig. (b). Under adiabatic conditions, this 
change in entropy is balanced by a decrease of tempera¬ 
ture, cf. Fig. The change of polarization and the cor¬ 
responding adiabatic temperature change are connected 
by Maxwell’s relation 

dT _ T (dP\ 

dE~~^ ■ 

The polarization increases by approximately 25 /rC/cm^ 
above Tq as the field of 100 kV/cm induces a large change 
of polarization in this temperature range, see Fig. 

In contrast, the change of polarization by an electrical 
field along [001] is one order of magnitude smaller within 
the ferroelectric phases at lower temperatures. Thus, 
the ECE due to the field induced ferroelectric transition 
above Tc of the heating-up simulation is largest and a 
maximal adiabatic temperature change of approximately 
3.7 K is obtained, see Fig. It has to be noted that 
the maximal adiabatic temperature change increases to 
7 K if pressure corrections are taken into account .1^ As 
discussed in Ref.[5|the ECE always increases with increas¬ 
ing T, cf. Eqn. ([^, and thus the increase of Tc by about 
150 K is the most important source of this modification. 

So far, all simulations have assumed an ideal bulk ma¬ 
terial. Local reduction of polarization will occur because 
of defects, such as oxygen vacancies, grain boundaries 
and surfaces. This setup is modeled for the limiting 
case of fully non-polar defects. These ’’inactive” dipoles 
reduce the overall dipole-dipole interaction in the sys¬ 
tem. Thus, the energy gain for a parallel alignment of 
all dipoles, i.e., a ferroelectric phase with large sponta¬ 
neous polarization, is systematically reduced. As the free 


energy (F) 


F{T, E) = Fo{T, E) - TS{T, E) 


( 3 ) 


rather than the energy (Fq), determines the equilibrium 
state of the system at finite temperatures, the transition 
temperature is shifted to lower temperatures, see Figs. 
and , since the paraelectric phase has a larger entropy 
(S). Also, the defects systematically reduce tetragonal 
strain (cf. Fig. 11), strain-dipole coupling energy, and the 


overall polarization compare lines with stars and squares 
in Fig. As the imposed defects may act as nucleation 
sites, they systematically reduce the thermal hysteresis, 
see Fig.|^ The steep change of the polarization at Tq and 
the first order character of the transition are conserved. 

As Tc and {dP/dT)E are reduced by the non-polar 
defects, the overall adiabatic temperature change (see 
Fig.[5| is systematically reduced with the number of de¬ 
fects, cf. Eqn. (i). For example, for AE = 100 kV/cm, 
5% non-polar defects reduce the ECE by a factor of two 
and shift the peak of the caloric response by approxi¬ 
mately 100 K to lower temperatures. The shape of the 
AT(T) peak is not considerably modified by the defects. 
Since especially no broadening of the curve is visible, the 
temperature range with finite temperature changes is lin¬ 
early reduced with the number of defects. 

In summary, non-polar defects reduce the ferroelectric 
phase transition temperature and polarization of BTO 
systematically. In turn also the maximal response of the 
ECE of BTO is reduced. However, even 5% defects do not 
suppress the caloric response and are thus no obstacle for 
the use of the material in ECE applications. In addition, 
the operation range and the temperature giving maximal 
adiabatic response may be lowered by non-polar defects. 
One has to note, that the present model overestimates 
the effect of local defects as ’’real” defects most likely 
correspond to a local reduction of the polarization rather 
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than a complete freezing and thus a smaller shift and 
reduction of the adiabatic response. 


4. POLAR COMPLEXES 

In the present section we discuss the influence of stable 
local defect dipoles. For example, this setup corresponds 
to dopant-oxygen vacancy complexes in transition metal 
doped perovskites. It has been reportecP^l^i^ that such 
complexes form immediately during the synthesis pro¬ 
cess, that they have much slower relaxation times than 
the free dipoles related to the ferroelectric soft mode, and 
that different relative orientations towards the sponta¬ 
neous polarization of the host ferroelectric material are 
possible. In the present study, we assume fixed defect 
dipoles, which is a rather realistic approximation at low 
and ambient temperatures during the time span of a typ¬ 
ical ECE measurement, cf. Ref. m 

Dipoles with a charging of 2 |e| and distance of 1.91 A 
have been found for fully equilibrated O vacancy-Cu com¬ 
plexes in PbTiOs.^ This would correspond to a local 
polarization of approximately 106 /xC/cm^ if one defect 
dipole per BTO unit volume is assumed. Fixed local soft 
mode amplitudes up to 0.4 A are taken into account, 
corresponding to a locally enhanced spontaneous polar¬ 
izations of 100 /iC/cm^ at the maximum. 


4.1. Unipolar cycling 

For ferroelectric perovskites, it has been shown that 
it is energetically most favorable if the defect dipoles 
order with the overall polarization of the ferroelectric 
phase.^i^^ Equilibration to this ground state (aging) can 
be achieved, e.g., by cooling in a strong external field. If 
a unipolar external electrical field is used in the ECE 
measurement, one may assume fixed defect dipoles par¬ 
allel to the direction of the field for the whole measuring 
time (switching on and off the field). 

The influence of 1% defect dipoles with strengths corre¬ 
sponding to soft mode amplitudes of 0.2 A (50 ^C/cm^) 
and 0.4 A (100 /xC/cm^) on the phase diagram with¬ 
out and with a parallel external field of 100 kV/cm are 
illustrated in Fig. With increasing strength or num¬ 
ber of defect dipoles the transition temperature is shifted 
to higher temperatures. This finding is important with 
respect to the interpretation of experimental results on 
doped BTO there either an increase, no modification, or 
even a decrease of Tc have been discussed. Our results 
show that fully equilibrated defects in a crystalline sam¬ 
ple may increase Tq. In contrast to this, defects may 
reduce Tc by extrinsic effects, such as the reduction of 
the grain size in experimental samples. 

Already for 1% defects with a soft mode amplitude of 
0.1 A the polarization profile is smooth and the thermal 
hysteresis is smaller than the accuracy of the simulation. 
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FIG. 6: (Color online) Polarization of BaTiOs with 1% po¬ 
lar defects aligned parallel to the external field of 100 kV/cm 
and after switching off the field. Black squares: no defects; 
blue open circles: 50/xC/cm^; red circles 100 /xC/cm^ (b) cor¬ 
responding approximate adiabatic temperature change. 


This soft mode amplitude corresponds to a local polar¬ 
ization of 25 /xC/cm^ which is even slightly smaller than 
the free dipoles of the host material. For an amplitude 
of 0.2 A, we find a rather smooth change of polarization 
with temperature, a vanishing thermal hysteresis, and a 
large polarization above the critical temperature of the 
ideal system. Thus, the system shows the characteristic 
features of a field induced polarization above the criti¬ 
cal field strength rather than the first order transition of 
the ideal system, see Fig. (a). The fixed local defect 
dipoles create an additional internal electrical field and 
thus the overall polarization and the transition tempera¬ 
ture increase with the strength of the defect dipoles and 
their density. 

Similarly, an internal bias field has been found exper¬ 
imentally in aged doped BTO.I^ There, the hysteresis 
is shifted along the field axis and the coercive field in¬ 
creases, if the defect dipoles in a doped sample have fully 
ordered. 

It has to be noted that the shift of the transition tem¬ 
perature by 100 K and the magnitude of the internal 
dipole field overestimate experimental results, e.g., in 
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Ref. [321 On the one hand, the used model for defects 
neglects relaxation effects and domains. One the other 
hand, it is not ensured that all dopants form fully relaxed 
dipole complexes in experiment, which anyway forbids a 
quantitative comparison of the bias field. In addition, 
we have assumed a maximal defect dipole amplitude of 
0.4 A corresponding to a rather large local polarization 
of 100 /rC/cm^ which may exceed the dipole moment of 
the real doped material. However, the same qualitative 
trends are found for a soft mode amplitude of 0.1 A. Here, 
the shift of Tc and thus the adiabatic response is only 
10 K. In summary, polar defects induce an internal elec¬ 
trical field which increases Tc and the polarization while 
the first order character of the ferroelectric transitions 
weakens. 

For the determination of the ECE the free dipoles are 
first equilibrated in an external field which is removed 
afterwards, see Secs. and During this field removal, 
the free dipoles relax between the two equilibrium phases 
found with and without external field as illustrated in 
Fig. [6l Compared to the ideal system without defects, 
the change of polarization between both states is reduced 
by the internal field induced by the defects. Although the 
system cools down under field removal, the fixed dipoles 
reduce the caloric response. In addition, the smoothing 
of the polarization profile reduces AT, cf. Eqn.[^ possible 
contributions by the latent heat of the first order tran¬ 
sition are lost as the system is beyond its critical point, 
and the number of free dipoles is slightly reduced by the 
defects. 

The caloric response of the ideal system approximately 
halves if 1% defects with a dipole moment correspond¬ 
ing to soft mode amplitudes between 0.1 and 0.4 A are 
imposed. We want to note that qualitatively the same 
caloric response is found by the indirect approach, i.e. by 
the integration of different equilibrium P(T,E) curves, 
see Sec. In the simulations, the maximal adiabatic 
temperature change and the shape of the AT(T) curve 
are basically the same for different dipole strengths (0.1- 
0.4 A) as soon as the internal field exceeds the critical 
field strength. The reduction of the peak maximum is 
accompanied by a slight broadening of the AT(T) curve. 
In turn a finite ECE if found for larger temperature range 
compared to the pure system. As the polarization pro¬ 
file is shifted to higher temperatures proportional to the 
strength of the internal dipole field, one could addition¬ 
ally adjust the optimal temperature of the caloric re¬ 
sponse by using proper dopants. Complementary to the 
dipole strength, the density of defect dipoles modifies the 
induced internal field. Thus, if the number of defects in¬ 
creases, the peak of the adiabatic response for an exter¬ 
nal field of 100 kV/cm is further reduced, broadened and 
shifted to higher temperatures. 

In summary, unipolar dipole defects pointing along the 
direction of the external field are no obstacle for a large 
ECE. Within the present model simulations an approx¬ 
imate ECE of 2 K can be found for an external field of 
100 kV/cm for 1% defects. In addition the AT(T) curves 
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FIG. 7: (Color online) Influence of antiparallel defects on 
the polarization of BaTiOs for different defect densities (4%: 
black curve at the bottom, 1%; otherwise), strengths (local 
soft mode amplitude of 0.4 A/0.2 A: stars/squares) and dif¬ 
ferent external fields (500 kV/cm: magenta; 50 kV/cm: cyan; 
100 kV/cm otherwise) for cooling-down simulations. For 1% 
defects with a strength of 0.4 A also results for heating-up 
(red) and initalized from scratch (black) as well as cooling- 
down simulation for parallel defects (green, with label P) are 
given. 

and thus the operation range of the ECE broadens, which 
is beneficial for operation. Only the narrow peak of ap¬ 
proximately 10 K width directly at Tc is reduced as the 
internal bias field is beyond the critical field. Although 
the used simple model may quantitatively overestimate 
the shift of the caloric response with temperature, the use 
of defect dipoles is clearly no obstacle and may even be 
used to tune the operation range of a ECE device under 
unipolar cycling. 

4.2. Antiparallel defects 

Depending on measurement protocol and history of the 
sample, one can also think of metastable defect dipoles 
pointing antiparallel to the overall po lariza tion for a cer¬ 
tain number of unipolar field cycles.li^E^ Such a setup 
may have exciting consequences for the phase diagram, 
see Figs.[^and[8| and the caloric response of the material, 
see Fig. 

We start our discussion with the influence of compet¬ 
ing internal (Tint) and external (Text) fields on the phase 
diagram of BTO. In the limit of an external field strongly 
exceeding the internal one, the free dipoles order along 
the external field direction (state 1), see Figs. and 
Tab. 1^ In this field range, the internal field acts as a 
small disturbance of the ferroelectric state, only. Small 
chains of free dipoles several lattice constants below and 
above the defect dipoles are coupled to the internal field 
resulting in a slight reduction of P and the tetragonal 
strain compared to the defect free material within an ex- 
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ternal field. With decreasing ratio between external and 
internal field, polarization and transition temperature of 
the ferroelectric state 1 are systematically reduced, see 
Fig.|8] 

In the opposite limit without external field, the inter¬ 
nal field stabilizes the ferroelectric phase with polariza¬ 
tion parallel to the defects (state 2), as discussed in the 
last section. With an increasing strength of the external 
field, this state (2) is destabilized, i.e. the polarization 
of the ferroelectric phase and its transition temperature 
systematically reduced, see Figs. and For compa¬ 
rable internal and external field strengths, the stability 
of the ferroelectric states 1 and 2 depends also on the 
history of the sample and the simulation protocol, see 
the example of a soft mode amplitude of 0.2 A in Fig. 
First of all, the competing fields open a thermal hystere¬ 
sis between cooling-down and heating-up simulations for 
state 2. Although, state 2 has the lowest energy and is 
found in cooling-down simulations up to an external field 
strength of 70 kV/cm, state 1 is at least a metastable 
state for E^xt A 30 kV/cm and is found for heating-up 
simulations or after field-removal. 


Below the transition temperature Tc, internal and ex¬ 
ternal field occasionally compensate each other, result¬ 
ing in a state without overall polarization (state 3). For 
example, for a soft mode amplitude of 0.4 A and an 
external field of 100 kV/cm a polarization of approxi¬ 
mately 2 fiC/cm^ is observed above 325 K/350 K (for 
cooling-down/ heating-up), which is twice the polariza¬ 
tion of pure frozen defect dipoles, see Fig. A snapshot 
of the local dipole configuration at 400 K reveals that 
the distribution of local soft mode amplitudes shows two 
maxima at approximately 17 and -13 ^C/cm^. Parts of 
the free dipoles align with external, while a similar num¬ 
ber of dipoles near the defects align with internal defect 
field. 

It has to be noted that all three states, their differ¬ 
ent temperature dependency and their dependency on 
the sample history can be reproduced with different ran¬ 
dom defect distributions, cell sizes, and defect strengths, 
see Tab. |l] In addition, the use of an external pressure 
does not modify the phase diagram despite the increase 
of transition temperature and polarization which as also 
been discussed for ideal BTO, cf. Figs. ^ and Tab. 

In summary, the stability region of three phases 
with pronounced differences in dipole ordering and 
polarization, can be adjusted by the relative strength of 
the external held and the internal held given by density 
and magnitude of the defect dipoles. In addition, 
diherent metastable states can be found depending on 
the simulation protocol and the history of the sample 
for similar strength of internal and external held. 


The rich phase diagram found for anti-parallel defects 
has exciting consequences on the ECE which will be dis¬ 
cussed in the following. Eirst of all, we want to highlight 
the dependency of the obtained phases on the history of 
the sample as discussed in the last paragraph. Because of 
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FIG. 8: (Color online) Polarization of BTO for 1% defects 
with a soft mode amplitnde of 0.2 A nnder cooling down (solid 
blue lines) and heating up (dashed red lines) for different field 
strength. In contrast to Fig. [^a temperature dependent ex¬ 
pansive pressure has been applied. A 16x16x16 supercell only 
has been used. 


TABLE I: Exemplary phase sequences found for different rel¬ 
ative strength of internal (JFint as given by dipole strength 
and concentration) and external fields (Aext) found for defect 
dipoles antiparallel to the external held in cooling-down and 


heating-up 

simulations. 




held 

phase sequence 

(a) 


cool: 

P = 0 => P||Pext 


heat: 

p||£’ext =4> P = 0 

e.g., 

1 % 0.4 A 500 kV/cm, 

1% 0.2 100 kV/cm 

(b) 

Eext Eint 

cool: 

P = 0^P||Pint 

heat: 

P||£ 

'ext P||Pint =4> P = 

e.g., 

1% 0.4 A 100 kV/cm, 

1% 0.2 A 70 kV/cm 

(c) : 


cool: 

P = 0 => P||Pint 


heat: 

P||Pint => P = 0 

e.g., 

1% 0.4 A 50 kV/cm, 4% 0.4 A 100 kV/cm , 

1% ( 

1.2 A 30 kV/cm 




this, also the caloric response depends crucial on the sam¬ 
ple history, the equilibration and the direction of the field 
change. As a consequence of this, the indirect method, 
i.e. the integration of one general set of equilibrium po¬ 
larization curves cannot be used similarly to the findings 
for relaxors.l^ 

We start our discussion with the caloric response un¬ 
der field-removal. First the system is equilibrated at each 
temperature within an external held. Approximately, the 
obtained phase corresponds to the heating-up phase di¬ 
agram, see Fig. If the external held is removed, it is 
energetically most favorable for the free dipoles to align 
parallel to the defects, i.e. if the state discussed in the last 
section for parallel dipoles is reached. However, with de¬ 
creasing temperature the relaxation slows down and thus 

























FIG. 9: (Color online) Approximate caloric response of BTO 
with 1% dipole defects if an external field of 100 kV/cm is 
removed. Black: no defects; cyan: parallel defects (u); Red: 
50:50 parallel and antiparallel defects (u/d); blue: antiparallel 
defects (d) with soft mode amplitudes of 0.4 A (solid lines) 
and 0.2 A (dashed lines); For the 0.4 A amplitude also the 
caloric response after field cooling (FC) is shown. 



T(K) 


FIG. 10: (Color online) Polarization for 1% polar defects 
with equal up/down distribution under different simulation 
protocols (cooling: blue, heating and after field removal: red) 
and within 100 kV/cm (black). For comparison P(T) under 
cooling is given for the ideal system (dotted lines). 

*: After field removal. 


the system with polarization along the field direction, 
may be stuck in this metastable state for typical simula¬ 
tion times. For example, after 500 ps the system relaxes 
to the energetic ground state at 210 K for defects with 
0.4 A amplitude after the external field of 100 kV/cm has 
been removed. 

For the fully relaxed system, field removal changes the 
polarization considerably. Especially, for initial states 1 
and 3, a large change of P under field removal induces a 
large temperature change, see Fig. For state 2, with 
a low coupling with Eext, the change of P and thus the 
adiabatic temperature change is reduced albeit a finite 
caloric response also appears in this temperature region. 

In contrast to the defect configurations discussed in 
Secs. and 4.1 the polarization is larger without field. 
Thus, an inverse ECE is found as the system heats up 
under field removal. The effect is rather large and espe¬ 
cially for the initial state 3 the magnitude of the caloric 
response is exceeding the conventional ECE found for 
pure BTO at the same strength of the external field, see 
Fig.ii It should be noted that the peak of the ECE 
around 350 K, related to state 3, does neither depend on 
the ramping rate, see discussion in Sec.[^ nor on the his¬ 
tory of the sample and is found for field-on and field-off 
simulations, see Eig. 13 (b), as the compensation between 
external and internal field is found uniquely between 
cooling-down and heating up simulations, cf. Figs. and 


Note, that for a dipole strength of 0.4 A also an inverse 
ECE of similar magnitude is found below 300 K, if an ex¬ 
ternal field of 100 kV/cm is removed from the initial state 
1. Similar AT values can even be found below 250 K for 
extended relaxation time (not shown in Fig. How¬ 


ever, for this combination of internal and external field, 
the metastable state 1 only appears for field-heated or un¬ 
poled samples and thus the ECE is not reversible and can 
only be found under field removal, cf. Fig. 13 (b). After 
the system has reached the stable state 2 with polariza¬ 
tion parallel to the defect dipoles in one of the first field 
cycles the caloric response is reduced below 1 K in the 
following ones, see blue stars in Fig. The temperature 
range with a large inverse ECE might be extended if the 
internal field strength is reduced, e.g., if a soft mode am¬ 
plitude of 0.2 A is assumed, see Fig.|^ In this case, states 
3 and 1 are stable within the external field of 100 kV/cm 
down to approximately 160 K, i.e. down to the second 
ferroelectric transition temperature between the tetrago¬ 
nal and orthorhombic ferroelectric phases. Uniquely, the 
same states are reached for cooling-down and heating- 
up simulations and thus a reversible effect has to be ex¬ 
pected. After field removal, the free dipoles align with 
the internal field induced by the defects if a sufficient re¬ 
laxation time is taken into account, cf. discussion for the 
soft mode amplitude of 0.4 A. As the polarization in the 
initial states is smaller than in the final state, an inverse 
ECE with similar magnitude as the conventional ECE of 
the ideal system is found in a broad temperature range 
below 350 K, see Fig. Even below 230 K a compara¬ 
ble caloric response has to be expected for an extended 
equilibration time after field removal. 


In summary, we find an inverse ECE in doped BTO if 
the external field and the internal field given by the defect 
dipoles are antiparallel. The relative strength between 
external field and internal field can be used in order to 
stabilize one or the other phase within the external field, 
thus adjusting the magnitude and even the sign of the 
ECE. Especially for similar magnitudes of external and 
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T(K) 

FIG. 11: (Color online) Illustration of the temperature de¬ 
pendent lattice constants along and perpendicular to the 
tetragonal axis for the ideal system (black), 1% non-polar 
defects (blue crosses), and 1% polar defects with 0.4 A am¬ 
plitude (dashed lines). Empty squares: no external field; 
Empty circles: 50:50 distribution of parallel and antiparal¬ 
lel dipoles; Filled circles: dipoles antiparallel to an external 
field of 100 kV/cm. Solid lines: cooling-down and dashed 
lines: after field removal. 


internal field (state 3) a reversible inverse ECE is found 
over a broad temperature range. The magnitude of this 
temperature change even exceeds the conventional ECE 
we find without defects. 

We note that the appearance of the large inverse ECE 
in doped ferroelectrics has been confirmed by simulations 
employing a Landau-type potential during the referee 
process of this paper.l^ 


under cooling for the configuration shown in Fig. 10 


as 


180° domains with up and downwards polarization are 
formed. However, the local soft mode, as well as the 
tetragonal ratio is very similar between cooling-down and 


heating-up, see Fig. 11 


One important mechanism for the shift of Tq can be 
found in the coupling between local strain and local soft 
mode. Thus, a tetragonal distortion is induced near 
the polar defects i.e. the lattice constant along the local 
dipole direction increases compared to the ideal system 
for all temperatures, see Fig. El This tetragonal distor¬ 
tion stabilizes the ferroelectric phase, the energy barrier 
for the ferroelectric transition is reduced, and Tc as well 
as the saturation polarization increase. If the in-plane 
lattice constants are clamped to their bulk value, the de¬ 
fect induced strain is reduced and also the shift of Tc to 
higher temperatures does not show up anymore in simu¬ 
lations of the material with frozen defects. 

As no net internal field builds up for the equal dis¬ 
tribution of both kind of defects, the polarization shows 
a steep jump at Tc for heating-up simulations or the 
ECE protocol. Analog to the ideal system, an exter¬ 
nal held induces a spontaneous polarization, resulting in 
large adiabatic temperature change above the transition 
temperature. However, the dipoles slightly reduce the 
adiabatic response at its peak maximum. This effect is 
similar to 1% non-polar dipoles but a slightly broader 
peak is found. 

In summary, for all relative orientations between lo¬ 
cal defects and external held, taken into account in this 
investigation, we hnd caloric responses which are compa¬ 
rable to the response of ideal BTO. Especially, the dipole 
distribution under an cycling external held is no obstacle 
for the ECE. The operation range is even broadened in 
comparison to the ideal system and it is slightly shifted 
to higher temperatures. 


4.3. Bipolar cycling 

Most device concepts so far are based on the bipolar 
cycling of an electrical held. For such a setup “de-aging” 
of the sample is present, i.e. the internal bias held dis¬ 
cussed for fully polarized samples with defects, vanishes 
with time. It has been discussed in literature that the 
slow defect dipoles cannot follow the electrical held and 
are thus exposed to a surrounding polarization pointing 
half times up and down. As consequence, also half the 
defect dipoles point in either held direction.!^ 

In the following, we assume hxed defect dipoles point¬ 
ing randomly either parallel or antiparallel to the exter¬ 
nal held with a 50:50 distribution. For 1% defects corre¬ 
sponding to a local soft mode of 0.4 A the transition tem¬ 
perature under cooling-down increases by approximately 
30 K compared to the ideal system, see Fig. [T0| 

Furthermore, the local defects act as kind of precursor 
and thus the thermal hysteresis between cooling-down 
and heating-up simulations is reduced to at most 5 K. It 
has to be noted that the overall polarization is reduced 


4.4. Disordered dipoles 

For the as-prepared cubic sample, the alignment of de¬ 
fect dipoles along all cubic lattice direction has the same 
energy and probability. In the following we discuss the 
inhuence on the phase diagram and the ECE by defect 
dipoles perpendicular to the applied held. 

The inhuence of a uniform distribution of 1% dipoles 
corresponding to a soft mode amplitude of 0.2 A along 
all cubic directions on the phase diagram is illustrated 
in Fig. |12| Heating-up and cooling-down simulations 
including a temperature dependent pressure correction 
have been performed. The jump of the polarization at Tc 
is steep and a thermal hysteresis is present, i.e. the hrst 
order character of the EEL transition is conserved. Here, 
only one sixth of the dipoles point parallel or antiparallel 
to the external held. This results in rather small internal 
helds along a certain direction. In addition, the dipole 
helds should in average cancel each other. In contrast to 
the example of dipoles pointing along and antiparallel to 
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TABLE II: Summary of the influence of different kind of defects on the phase diagram and caloric response as obtained within 
our ab initio based molecular dynamics simulations. 


Impurity 

non-polar defect or 
locally reduced polarization 

Tc 

reduced 

ECE 

reduced 

no broadening 

thermal hysteresis 

reduced 

polar defects 
parallel to external field 

increase 

reduced 

broadening of AT(T) 

vanishes 

50:50 distribution of 

parallel and antiparallel defects 

increase 

reduced 

broadening of AT(T) 

vanishes 

antiparallel defects 
in unipolar field 

increase 

inverse 

broad range 

different metastable 

states 

randomly ordered 
dipoles 

reduced 

reduced 

reduced 



T(K) 


FIG. 12: (Color online) Cooling-down and heating-up simu¬ 
lation including a temperature dependent pressure correction 
for a 16x16x16 simulation cell. 1% defects corresponding to 
a soft mode amplitude of 0.2 A have been taken into account 
for dipoles pointing randomly along the six cubic crystal di¬ 
rections (blue diamonds) and and non-polar defects (black). 
The polarization profiles without external held (solid lines) 
and in an external held of 100 kV/cm along z (dashed lines) 
are opposed. 


the overall polarization, no direction for the tetragonal 
strain is favored by the defect dipoles and thus Tc is not 
shifted to higher temperatures. Instead, the magnitudes 
of dP/dTE and Tc are reduced by the doping similar to 
1% non-polar defects and thus a slightly reduced ECE 
has to be expected. 

In summary, the equal distribution of defect dipoles 
along all cubic lattice directions in an as-prepared cu¬ 
bic sample may slightly reduce the ECE. However after 
several ECE cycles, the defect dipoles would most likely 
equilibrate with the direction of the applied fielcP^^hl! and 
the ECE as discussed in Sec. |4] would recover. 


5. CONCLUSIONS AND OUTLOOK 

We have discussed the influence of different defects on 
the ferroelectric and electrocaloric properties of BTO in 
the framework of ab initio based molecular dynamics sim¬ 
ulations. The main influences of different kinds of defects 
are summarized in Tab. ttn We find that defects with 
locally reduced polarization tend to reduce Tq and the 
adiabatic response. The same trend is found for strong 
defect dipoles with slow relaxation, but without a pre¬ 
ferred alignment. However, in both cases the reduction 
of the ECE is only small for the tested defect concentra¬ 
tions and thus does not contradict a use of BTO in ECE 
cooling devices. 

In addition, we have discussed the possible benefit of 
dopants for the ECE. Eor this purpose we have simu¬ 
lated the influence of fixed defect dipoles on the ferroelec¬ 
tric phase diagram and the ECE systematically. Dipoles 
pointing along the tetragonal axis of BTO tend to in¬ 
crease Tc which could partly be related to the strong 
dipole-strain coupling. In addition, parallel dipoles in¬ 
duce an internal electrical field which results in an in¬ 
creased polarization and reduces the first order character 
of the ferroelectric transition. 

If the defect dipoles, or at last half of them, are point¬ 
ing parallel to the external field, the caloric response is 
slightly reduced and the AT curve is shifted to higher 
temperatures and broadens. This kind of defect dipoles 
may be beneficial for applications based on ferroelectric 
materials as they open up a possibility to broaden the 
operation range. Excitingly, the doping of BTO samples 
may even induce different stable and metastable ferro¬ 
electric states and an inverse ECE under certain mea¬ 
surement protocols. Our investigation has shown that 
defect dipoles pointing antiparallel to the external field 
form an internal electrical field which couples to the free 
dipoles with a different temperature dependency than 
the external field. The relative strength between internal 
and external field at a certain temperature may be ad¬ 
justed, e.g., by the number and strength of the imposed 
defect dipoles. As a consequence three different ferroelec- 
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trie states with polarization along external field, internal 
field, or nearly compensated polarization can be seen in 
our simulations. 

If the field is removed, the free dipoles align again 
with the defect dipoles resulting in an increasing polar¬ 
ization and thus a heating of the sample for all these 
states. Especially the high-temperature state which we 
found for different combinations of internal and external 
field is promising for the modification of the ECE: In our 
simulations the polarization within the external field is 
quenched in a broad temperature interval below T^, re¬ 
sulting in a large change of the polarization and a large 
caloric response under field removal. 

Euture work will concentrate on improved modeling 
of the influence of defects, e.g. by using slowly relaxing 
but dynamic defect dipoles and by taking the coupling 
between external field and defect dipole strength into ac¬ 
count. In addition, the combination of defects with alloy¬ 
ing is promising with respect to adjusting the operation 
temperature for possible applications. 
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6 . APPENDIX 

In the following, technical aspects of the simulation 
will be discussed in detail, a comparison of the direct 
and indirect method will be reviewed from literature,^ 
and additional convergency tests will be presented. In 
the present paper we calculate the caloric response of 
BTO directly, see Sec. For this purpose we equilibrate 
the system at a given temperature (Ti) and external field 
£i, ramp the field to its final value £ 2 , and monitor the 
change of the kinetic energy (T 2 ). By this approach we 
can get a good qualitative description of the change of 
the ECE with different kind of defects as we take the 


most relevant degrees of freedom into acc ount and i nduce 
similar errors without all our simulations! ^ * ^^ * ^^ *^ 

Alternatively, the thermodynamic Maxwell relations 
(cf. Eqn. (§) can be used to calculate the caloric re¬ 
sponse by integration of the pyroelectric response at dif¬ 
ferent values of the external field (indirect method) 


AT = - 



T 



( 4 ) 


Here, Cp^£ is the specific heat at constant pressure and 
field, and the external field is varied form £i to £ 2 - This 
indirect approach breaks down at the first order phase 
transition as dP/dT is ill-defined, the specific heat di¬ 
verges, and possible contributions stemming from the la¬ 
tent heat are not accounted for.E^ In addition, the indi¬ 
rect method fails if the caloric response depends on the 
simulation protocol, the history of the sample, or the 
direction of the field change (on-off), e.g., in relaxors.^^ 
However, the indirect approach has the important advan¬ 
tage that only equilibrium states are taken into account 
and it can be used in order to cross check the results 
obtained via the direct approach, which include two im¬ 
portant approximations. 

First, in order to compare our calculated tempera¬ 
ture changes with experimental results, one has to keep 
in mind that we treat only three degrees of freedom 
as dynamical variables during our MD simulations, see 
Sec. During the time evolution the algorithm thus 
takes the system from the state (Pi,Ti,£’i), to the new 
state {P 2 ,T 2 ,£ 2 ) while the state for the actual number of 
degrees of freedom would be {P'i,T^,£ 2 )- The leading er¬ 
ror in the adiabatic temperature change can be corrected 
if we rescale our results with a factor of 15/3, i.e. by 
rescaling the number of the degrees of freedom. Due to 
the uncorrected error in the final polarization (P 3 ), an 
error smaller than 1 K is i nduce d in the calculated adi¬ 
abatic temperature changewhich does not modify 
the trends of the caloric response for different kinds of 
defects. 

Second, experimental ramping rates are not accessible 
in our simulations due to the computational effort. For 
pure BaTiOa, extensive convergency tests with respect to 
the ramping rate have been performed in Ref. 1191 show¬ 
ing full convergency for a ramping rate of 0 . 002 kV/cm/fs. 
In addition, an errorbar of about 0.3 K has been found 
for the ramping rate of 0.05 kV/cm/fs which is used in 
the present investigation. Indeed, the adiabatic temper¬ 
ature change of BTO as obtained by the direct method 
with these ramping rates and the adiabatic temperature 
change as obtained by the indirect method coincide as 
soon as £i and £2 are beyond the critical field strength, 
cf. Ref. [ini Thus, the reliability and accuracy of the 
used direct method is fully supported by the compari¬ 
son with the indirect one. Furthermore, our calculated 
temperature change is in good qualitative agreement to 
experimental values, e.g. AT =0.9 K has been found for 
an applied field of 12 kV /enP^ while we yield AT =3.8 K 
for an applied field of 100 kV/cm. 
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FIG. 13: (a) Comparison of the adiabatic temperature change 
found by the direct (dotted lines) and indirect (solid lines) 
method for 1% parallel defects with a strength corresponding 
to a soft mode amplitude of 0.2 A. Black: Field range: 0— 
100 kV/cm; Red: 0—50kV/cm. Lines with symbols: ramping 
down the field with a rate of 0.05kV/cm/fs otherwise instan¬ 
taneous ramping has been used, (b) Corresponding entropy 
change found by the indirect method, error bar 


In order to test the influence of defects on the accu¬ 
racy of the direct approach, the rescaling of the degrees 
of freedom, and the fast ramping, we have performed a 
set of calculations for 1% polar defects parallel to the ex¬ 
ternal field with a strength corresponding to 0.2 A, see 
Fig. 13 (a). First, we have calculated the adiabatic re¬ 
sponse by ramping down a field from £i = 100 kV/cm 
to £2 = 0 kV/cm with a field rate of 0.05 kV/cm/fs. 
Second, we have switched off the field instantaneously. 
Third, the indirect method has been applied, using the 
high-temperature value of C'p_£=15 k^ =3.3 J/Kcm^ as 
quantum mechanical effects are neglected within our sim¬ 
ulations. We note that the indirect method can be ap¬ 
plied even for £2 =0 kV/cm as the system is beyond its 
critical point due to the strong internal field induced by 
defect dipoles, see Sec. |4.1[ 

Although, instantaneous switching may indeed bring 
the system out of equilibrium, see the detailed discus¬ 


FIG. 14: (Color online) Dependency of the adiabatic response 
on the ramping rate for 1% antiparallel defects with a strength 
corresponding to a soft mode amplitude of 0.4 A. The ramping 
rate is given in kV/cm/fs. Black and red: field removal, Blue: 
field on. 


sion in Ref. m the quantitative change of the caloric 
response compared to held ramping is already well be¬ 
low (below 0.25 K) the different trends of the adiabatic 
response discussed throughout this paper. In addition, 
there is a rather good qualitative and semi-quantitative 
agreement between direct and indirect approach. Small 
discrepancies can be understood by means of the temper¬ 
ature dependence of the specific heat. While the specific 
heat is approximately constant far from the phase transi¬ 
tions, there are pronounced peaks around Tc- Thus, es¬ 
pecially at the maximum of the AT peak and at slightly 
lower temperatures, AT as obtained with the indirect 
method is too large due to the approximation used for 
Cp^e- For example, if the specific heat calculated for the 
ideal system for an intermediate field strength, such as 
75 kV/cm is used, the peak maxima between direct and 
indirect approach coincide, see the detailed discussion in 

Ref.m 

In summary, both direct and indirect method yield the 
same results, despite the different approximates made, 
underlining the good qualitative estimation of the tem¬ 
perature change found in our simulations and the small 
influence of possible errors discussed above. 

In addition to the adiabatic temperature change (AT), 
the isothermal entropy change (AS) can be used to quan¬ 
tify the ECE. In the indirect approach this quantity is 
given by the Maxwell relation as 




d£ 


( 5 ) 


Figure 13 (b) exemplary shows the isothermal entropy 
change if an external held is removed from the sam¬ 
ple with 1% parallel defects. For £i = 100 kV/cm 
and £2 = 0 kV/cm we find a maximal increase in en¬ 
tropy of AS'=4.5 J-kg“^K“^. This value is comparable 
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to AS = 7 J-kg ^ found for ideal BTO and a field 
range from £i = 300 kV/cm to £2 = 75 kV/cmP^ It has 
to be noted that AS" and AT are approximately related 
to each other by 


AT - -^A5'E 

^P,€ 


( 6 ) 


Thus, the same trends of AS' and AT have to be expected 
for the various defect configurations discussed through¬ 
out this paper. 

For antiparallel defects, the use of the indirect methods 
may fail, as the phase diagram depends on the history of 
the sample, e.g. cooling- down or heating-up simulations, 
see Fig. Instead, the convergency of the adiabatic 
temperature change with respect to the ramping rate is 
illustrated for antiparallel defects in Fig. M 1% dipoles 
corresponding to 0.4 A have been used and an external 
field of 100 kV/cm has been removed or applied with 


two different ramping rates. Obviously, the adiabatic 
temperature change down to starting temperatures of 300 
K is sufficiently converged with respect to the ramping 
rate and a reversible temperature change is found for 
ramping the field on and off. Analogous, the adiabatic 
temperature change for 1% defect dipoles corresponding 
to 0.2 A as obtained with ramping rates of 0.5 or 0.001 
kV/cm/fs differ less than 0.1 K for temperatures above 
300 K. 


We note, that for the different metastable states, which 
occasionally appear at lower temperatures or different 
combinations of internal and external field strength, the 
direction and ramping rate of the field may play a role. 
For example, for 1% defect dipoles corresponding to 0.2 A 
AT as obtained with different ramping rates below 300 K 
differs by 0.5 K. We leave a closer inspection of this pa¬ 
rameter range for upcoming studies. 


^ X. Moya, S. Kar-Narayan, and N. D. Mathur, Nature 
Mater. 13, 439 (2014). 

^ M. Marathe and C. Ederer, App. Phys. Lett. 104, 212902 
(2014). 

® Y. Bai, G.-P. Zheng, K. Ding, L. Qiao, and S.-Q. Shi, J. 

Appl. Phys. 110, 094103 (2011). 

^ S. P. Beckman, L. F. Wan, J. A. Barr, and T. Nishimatsu, 
Mat. Letters. 89, 254 (2012). 

® M. Ozbolt and A. P. A. Kitanovski and, J.Tusek, Interna¬ 
tional Journal of Refrigeration 37, 16 (2014). 

® 1. Ponomareva and S. Lisenkov, Phys. Rev. Lett. 108, 
167604 (2012). 

^ J. F. Scott, Annu. Rev. Mater. Res. 41, 229 (2011). 

® S. Kar-Narayan and N. D. Mathur, J. Phys. D: Appl. Phys. 
43, 032002 (2010). 

® H.-J. Hagemann, J.Phys.C.:Solid State Phys: 11, 3333 
(1978). 

M. Maglione, R. Bohmer, A. Loidl, and U. T. Hochli, Phys. 
Rev. B 40, 11441 (1989). 

R. Maier, J. L. Cohn, J. J. Neumeier, and L. A. Bendersky, 
App. Phys. Lett. 78, 2536 (2001). 

H. Neumann and G. Arlt, Ferroelectrics 76, 303 (1987). 

C. H. Park and D. J. Chadi, Phys. Rev. B 57, R13961 
(1998). 

Y.-K. Choi, T. Hoshina, H. Takeda, and T. Tsurumi, Jpn. 
J. Appl. Phys 50, 031504 (2011). 

P. Erhart, P. Traskelin, and K. Albe, Phys. Rev. B 88, 
024107 (2013). 

P. Erhart, R.-A. Eichel, P. Traskelin, and K. Albe, Phys. 
Rev. B 76, 174116 (2007), URL http://link.aps.org/ 
doi/10.1103/PhysRevB.76.174116 

P. Erhart and K. Albe, J. Appl. Phys. 102, 084111 (2007). 

S. Lisenkov and 1. Ponomareva, Phys. Rev. B 80, 
140102(R) (2009). 

M. Marathe, A. Griinebohm, T. Nishimatsu, P. Entel, and 
C. Ederer, Phys. Rev. B 93, 054110 (2016). 

T. Nishimatsu, U. V. Waghmare, Y. Kawazoe, and D. Van¬ 
derbilt, Phys. Rev. B 78, 104104 (2008), URL http: 
//loto.sourceforge.net/feram/ 

R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 49, 


5828 (1994). 

W. Zhong, D. Vanderbilt, and K. M. Rabe, Phys. Rev. B 
52, 6301 (1995). 

T. Nishimatsu, M. Iwamoto, Y. Kawazoe, and U. V. Wagh¬ 
mare, Phys Rev B 82, 134106 (2010). 

S. D. Bond, B. J. Leimkuhler, and B. B. Laird, 

Journal of Computational Physics 151, 114 (1999), 

ISSN 0021-9991, URL http://www.sciencedirect.com/ 
science/article/pii/S002199919896171X 

T. Nishimatsu, J. A. Barr, and S. P. Beckman, J. Phys. 
Soc. Jpn. 82, 114605 (2013), URL http://jpsj .ipap.jp/ 
link?JPSJ/82/114605/ 

®® W. Zhong, D. Vanderbilt, and K. M. Rabe, Phys. Rev. 
Lett. 73, 1861 (1994). 

O. Dieguez, S. Tinte, A. Antons, C. Bungaro, J. B. Neaton, 
K. M. Rabe, and D. Vanderbilt, Phys. Rev. B 69, 212101 
(2004). 

S. Lisenkov and 1. Ponomareva, Phys. Rev. B 86, 104103 

( 2012 ). 

B.-K. Lai, 1. Ponomareva, 1. Kornev, L. Bellaiche, and 
G. Salamo, Appl. Phys. Lett. 91, 152909 (2007). 

J. Paul, T. Nishimatsu, Y. Kawazoe, and U. V. Waghmare, 
Phys. Rev. Lett. 99, 077601 (2007). 

A. Kumar and U. V. Waghmare, Phys. Rev. B 82, 054117 

( 2010 ). 

T. Kanata, T. Yoshikawa, and K. Kubota, Solid State 
Communications. 62, 765 (1987). 

G. Arlt and H. Neumann, Ferroelectrics 87, 109 (1988). 

S. Lu, B. Rozic, Q. M. Zhang, Z. Kutnjak, R. Pric, M. Lin, 

X. Li, and L. Gorny, Appl. Phys. Lett 97, 202901 (2010). 

Y. -B. Ma, K. Albe, and B.-X. Xu (2015), URL http:// 
arxiv.org/pdf/1507.05004vl.pdf 

X. Moya, E. Stern-Taulats, S. Crossley, D. Gonzalez- 
Alonso, S. Kar-Narayan, A. Planes, L. Manosa, and N. D. 
Mathur, Advanced materials 25, 1360 (2013). 

It has to be noted that instantaneous field removal and an 
explicit treatment of the acoustic degrees of freedoms have 
been used in connection with the pressure correction. 





